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When described through a plane- wave basis set, the inclusion of exact nonlocal exchange in 
hybrid functionals gives rise to a singularity, which slows down the convergence with the density of 
sampled k points in reciprocal space. In this work, we investigate to what extent the treatment of 
the singularity through the use of an auxiliary function is effective for fc-point samplings of limited 
density, in comparison to analogous calculations performed with semilocal density functionals. Our 
analysis applies for instance to calculations in which the Brillouin zone is sampled at the sole F-point, 
as often occurs in the study of surfaces, interfaces, and defects or in molecular dynamics simulations. 
In the adopted formulation, the treatment of the singularity results in the addition of a correction 
term to the total energy. The energy eigenvalue spectrum is affected by a downwards shift of the 
energy eigenvalues of the occupied states, while those of the unoccupied states remain unaffected. 
Analogous corrections also speed up the convergence of screened exchange interactions despite the 
absence of a proper singularity. Focusing first on neutral systems, both finite and extended, we show 
that the account of the singularity corrections bears convergence properties which are quantitatively 
similar to those observed with semilocal density functionals. We emphasize that this is not the 
case for uncorrected energies, particularly for elongated simulation cells for which qualitatively 
different trends are found. We then consider differences between total energies of systems differing 
by their charge state. For systems involving localized electron states, such as ionization potentials 
and electron affinities of molecular systems or charge transition levels of point defects, the proper 
account of the singularity correction yields convergence properties which are similar to those of 
neutral systems. In the case of extended systems, such energy differences provide an alternative 
way to determine the band edges, but are found to converge more slowly with simulation cell than 
in corresponding semilocal functionals because of the exchange selfinteraction associated to the extra 
charge. 

PACS numbers: 71.15.Mb 71.55.-i 



I. INTRODUCTION 

In the past decades, density functional theory-^ has 
become the mainstream technique for electronic struc- 
ture calculations of large molecules, clusters, liquids, and 
solids. In condensed matter applications, the most com- 
mon functional has initially been based on the local den- 
sity approximation (LDA)^ but generalized gradient ap- 
proximations have become increasingly popular in the 
last two decades £d£JLil A class of functionals that could 
potentially lead to higher accuracy include a fraction 
of exact nonlocal exchange in the exchange-correlation 
potentiali^i^ This class of functionals, referred to as hy- 
brid functionals, has become the standard electronic- 
structure approach in quantum-chemistry applications. 
For molecular systems, these functionals achieve a more 
accurate description not only of atomization energies^ 
but also of ionization potentials and electron affinities. 11 
However, several shortcomings still remain. For instance, 
the desired chemical accuracy is not always attained and 
no solution is offered for the treatment of Van-der-Waals 
interactions. Nevertheless, it appears quite clearly that 
the inclusion of exact exchange constitutes an improve- 
ment, which might be particularly useful in several cir- 
cumstances. 

When applied to semiconductors and insulators, hy- 



brid functionals provide a superior description to semilo- 
cal functionals. For instance, structural parameters are 
found to be closer to experimental values. 12 ' 13 Further- 
more, hybrid functionals give electronic band gaps which 
are systematically larger than those achieved with semilo- 
cal functionals, generally leading to a better agreement 
with experiment i 12 i 13 ' 14 In particular, the improvement 
achieved by hybrid functionals in which the Coulomb 
potential is screened is remarkable ; 12 : 13 and the ori- 
gin of this successful description can be rationalized. 15 
A more accurate description of the band gap is espe- 
cially important in certain classes of problems, such 
as the study of surface and interface state o 16 ' 17 and 
the determination of defect energy levels^ Indeed, a 
physically meaningful description of electronic states 
lying in the band gap can only be achieved when 
the calculated band gap approaches the experimental 

One 17 - 18 - 19 - 20 - 21 - 22 - 23 - 24 ' 25 - 26 -27.28,29,30,31,32.33,34 

For the treatment of condensed systems, hybrid func- 
tionals have been available for some time in codes 
combining Gaussian basis sets and periodic boundary 
conditions. 35 ' 36 However, it is expected that the treat- 
ment of the electronic structure through hybrid func- 
tionals carries potential to be much more widely used 
in an implementation based on plane-wave basis sets and 
pseudopotentials; 3 ? In this respect, the treatment of ex- 
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act exchange poses several new problems with respect to 
the use of standard semilocal functionals. 

First, the calculation of exact exchange entails a signif- 
icantly higher computational cost. To address this issue, 
efficient algorithms have been developed to optimize the 
scaling with respect to the number of plane waves . 38 ' 39 ! 40 
A further gain is achieved through optimal adaptation 
to new massively parallel computer platforms^ Despite 
these improvements, plane-wave based hybrid-functional 
calculations for a system of 1000 electrons still yield 
a computational cost exceeding that of semilocal ones 
by more than two orders of magnitude. 28 ' 42 Neverthe- 
less, such calculations are attracting increasing interest 
as larger computational resources become available. 

Second, the expression for exact exchange includes an 
integrable divergence,— which hinders its straightforward 
use within plane-wave formulations because of its slow 
convergence with the density of k points. Gygi and 
Baldereschi proposed a numerical treatment of the di- 
vergence based on the analytic integration of an auxil- 
iary function showing the same singularity 44 Such aux- 
iliary functions are nowadays available for arbitrary unit 
cells . 39 ' 45 ' 46 Alternative treatments consist in truncat- 
ing the Coulomb operator,— using a screened Coulomb 
potential)* - or transforming the Bloch functions in order 
to compute real-space Coulomb integrals. 49 

Third, the nonlocal exchange coupling between va- 
lence and core states also intervenes in the basic pseu- 
dopotential approximation. Such interactions can be ac- 
counted for within an all-electron scheme in which the 
pscudopotential approximation consists in freezing the 
wave functions of the core states^ The core-valence in- 
teractions due to exchange then lead to a modification 
of the nonlocal pseudopotential term. In consideration 
of the fact that current hybrid functionals generally only 
include a fraction of exact exchange (about 25%), pseu- 
dopotentials derived within semilocal formulations have 
often been transferred to hybrid schemes without any 
modificatio n 28 ' 29 ' 40 ' 50 or through the only use of nonlin- 
ear core corrections^ However, when core- valence inter- 
actions are sizable, these practices require particular care, 
and it has recently been pointed out that they might lead 
to inconsistencies^ 

In this work, we focus specifically on aspects associated 
with the treatment of the integrable singularity when 
calculating the exact exchange expresssion in electronic- 
structure schemes based on plane-wave basis sets. Fol- 
lowing Gygi and Baldereschi, 44 - we adopt a formulation in 
which the divergence is treated analytically. This scheme 
can be recast in such a way that the treatment of the di- 
vergence results in a correction term formally correspond- 
ing to the Fourier component of the exchange potential 
at vanishing wavevector in the Brillouin zone.— We first 
address the degree of convergence that can be achieved 
with a fc-point sampling of finite density for both the 
total energy and the energy eigenvalues, in comparison 
with a similar calculation that does not include exact 
exchange. This aspect is particularly important to vali- 



date calculations based on large simulations cells with the 
Brillouin zone sampled only at the T point, a configura- 
tion which is often used in surface, interface, and defect 
calculations or in ab initio molecular dynamics simula- 
tions. We also consider exchange interactions based on 
screened Coulomb potentials. While this case formally 
does not show a singularity, there are specific circum- 
stances in which the same convergence deficiency occurs 
as for the bare Coulomb potential. We then study the 
effect of the singularity correction on the total energy of 
charged systems. Differences between total energies of 
systems involving a different number of electrons are rel- 
evant for determining the electron affinity, the ionization 
potential, the band gap, and the charge transition levels 
of defects. We discuss the effect of the correction on these 
quantities for the cases of both finite and extended sys- 
tems. In particular, in the case of extended systems, we 
clarify the way the singularity correction affects the band 
edges obtained from the energy eigenvalues and those ob- 
tained through total-energy differences. 

This paper is organized as follows. In SecHH we review 
the treatment of the singularity of the exact exchange op- 
erator in the case of plane-wave basis sets following the 
scheme of Gygi and Baldereschi^ 4 . An extension of this 
method to the case of large supercells with sparse fc-point 
sampling is described. A generalized procedure for treat- 
ing the case of screened exchange is also given. In Sec. 
IIII1 we describe the computational setups used in this 
work. Section IIVI presents convergence tests for total 
energies, energy eigenvalues, and single-particle energy 
gaps of neutral systems including molecules and solids. 
Results obtained with T-point sampling are studied for 
increasing supercell size. In the case of extended sys- 
tems, a comparison is carried out with converged calcu- 
lations achieved with small unit cells and dense fc-point 
samplings. In Sees. IVl and IVI1 we study charged systems 
with localized and delocalized charge states, respectively. 
In particular, we focus on total energies, single-electron 
eigenvalues, electron addition energies, electron removal 
energies, and charge transition levels of defects. We study 
the convergence of these quantities as a function of the 
singularity correction and emphasize the difference be- 
tween localized and delocalized states. The conclusions 
are drawn in Sec. IVII1 

We note that a method related to hybrid function- 
als consists of including exact exchange in a local way 
through the use of an optimized effective potential. 52 ' 53 
While this scheme is not explicitly discussed in the fol- 
lowing, our considerations regarding the singularity cor- 
rection also apply in this case. 



II. TREATMENT OF THE SINGULARITY 

A. Exact exchange 

The exchange energy of a solid is a finite quantity if 
expressed per one unit cell 4 s - In a basis of plane waves, 
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the matrix element of the Fock exchange operator is given by^ 4 - 
I 
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where the sum over m runs over the occupied states and the integral is carried out over the Brillouin zone (BZ). 
Since the integrand diverges for q = k, special numerical care is required when replacing the integral with a sum 
over a finite number of fc-points. To treat the singularity, Gygi and Baldereschi used an auxiliary function, periodic 
in reciprocal space and showing the same ~ l/(k — q) 2 divergence as the integrand in Eq. (fT]). The integral of this 
function is then subtracted and added to the right hand side of Eq. JT]) . The subtracted term eliminates the divergence 
of the integrand and turns it into a smooth function of q, which can accurately be evaluated through a sampling of 
special fc-points. The singularity is effectively transferred to the added term and is taken care of through analytical 
integration. 44 ' 46 ! 47 

The method of Gygi and Baldereschi requires to be adapted in order to be applied to calculations with large 
supercells and sparse fc-point samplings. First, it is convenient to adopt the notation Q = q+G" and c mq (G + G") = 
c m (G + Q): 
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where the integral is over the whole of reciprocal space. For illustration, we focus in the following on a fc-point 
sampling based on the sole T point, i.e. k = T. Application of the procedure proposed by Gygi and Baldereschi 
transforms the right hand side of Eq. @ to 



c* m (G' + Q)c m (G + Q) 

Q 2 



c* m (G')c m (G)/(Q) 



-Lj2 c *m(G')c m (G) Jf(Q)dQ, 



(3) 



where the auxiliary function needs to chosen in such a way that /(Q) — > 1/Q 2 when Q — > 0. In addition, the function 
/ must be integrable in the whole of reciprocal space. The term in parentheses is now a smooth function of Q in 
which the divergence is cancelled. Therefore, it is justified to approximate the integral over Q in the first term in Eq. 
([3]) with a discrete sum via 
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E> 



(4) 



where f2 is the volume of the simulation cell. We here assume that the discretization in Q-space corresponds to the 
fc points for which the wave functions are determined. The second term in Eq. ^ is to be evaluated through an 
analytical integration. The final expression reads: 
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-L^c r *„(G') Cm (G) J f (Q) dQ 



In the case of T point sampling, the sum over Qi simply corresponds to a sum over the reciprocal lattice vectors G" 
from which G" = needs to be excluded. I 



We emphasize that the closer /(Q) is to 1/Q 2 , the 
smoother the integrand in the parentheses, and the more 
accurate the approximation of the integral by a sum of 
discrete terms. Hence, this condition should guide the 
choice of the optimal function /. We note that in the for- 
mulation of Gygi and Baldereschi, the auxiliary function 
is to a certain extent arbitrary, because ultimate conver- 
gence can always be achieved by increasing the density 



of the fc-point sampling. In the case of a fixed fc-point 
sampling, the converged result is achieved by consider- 
ing simulation cells of increasing size. However, in prac- 
tice, this limit is computationally more prohibitive. For 
instance, as a word of caution, we emphasize that for 
small band-gap materials very large supercells might be 
required in calculations relying only on the T-point even 
at the semilocal level. Generally, a satisfactory target for 
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calculations including exact exchange corresponds to the 
achievement of the same level of convergence as that at- 
tained for a semilocal functional under the same /c-point 
sampling conditions. For this reason, it is important 
to chose a function / which yields the best convergence 
properties for the chosen fc-point sampling. 

By reorganizing the order of the terms in Eq. ([5]), we 
obtain the following appealing form: 



G")c m (G + G")$(G"), (5) 
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£ c » ( G ' - 



,,G" 



where $ (G) represents a suitable generalization of the 
Fourier transform of the exchange potential and is given 
by: 



{1 4"7T 
X for G = 0, 
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This form is particularly convenient since it implies that 
only the G = term needs to be modified in practical im- 
plementations. We note that standard numerical imple- 
mentations do not experience any difficulty for treating 
the differences between large numbers that this reorga- 
nization of terms implies. 

A suitable form of the auxiliary function / is : 39 ' 45 



/(Q) = 



-lQ 2 



(8) 



For this particular choice of auxiliary function, the G=0 
term of the exchange potential becomes: 
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For illustration, we show in Fig. Q] the dependence of \ 
on 7 in the case of a cubic cell with a side of 20 bohr. As 
mentioned above, the approximation of replacing the in- 
tegral in Eq. with a discrete sum is the more accurate, 
the closer /(Q) is to 1/Q 2 . This leads to the following 
well-defined expression for x : 



x = x(o) 
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We now analyze the effect of the singularity correc- 
tion on eigenvalues and total energies. In our notation, 
the superscript "corr" indicates that the correction % nas 
been accounted for, whereas the superscript "uncorr" is 
used for quantities obtained with <I>(G=0) = in the 
definition (0. 
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FIG. 1: The function \ [Eq. (H2 vs the parameter 7 for a cubic 
simulation cell with a side of 20 bohr and F-point sampling. 
The optimal singularity correction is achieved for 7 = 0. 



The exact exchange term contributes to the eigenvalue 
e n of the state n through the following expression: 



Ae„= £ <(G)c„(G')(G 

G,G' 

where we used the matrix element ( G 
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V x G'y given in 

Eq. ([5]). The singularity correction affects the eigenvalue 
as follows: 

a corr a uncorr 

= -X £ £ < (G) c„ (G') c* m (G) c m (G') 

G,G' m 



"X£^m = I 
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where m runs over the occupied states and where we 
used the orthonormalization condition of the eigenstates. 
Thus, all the eigenvalues of occupied states shift down by 
X, whereas those of unoccupied states remain unaffected: 



X for n occupied, 
for n unoccupied. 



(13) 



In this derivation, we assumed that the ordering of the 
states is not affected by the application of the singularity 
correction. 

The exact exchange part of the total energy is given 

by 



E * = ?£ £ <(gk(g')(g 



F x 



G' 



(14) 



n G,G' 

where the contribution of the singularity correction is 
-f £ £<(GK(G')C (G)c m (G') = 

G.G' n.rn 

with N c \ corresponding to the number of electrons in the 
supercell. Thus: 



pcorr puncorr 

-^tot — ^tot 



(16) 



i.e. a correction of — x/2 for each electron. 
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In the limit 7 — ► [cf. Eq. (|10p]. the correction — x/2 
corresponds to the electrostatic energy of a point charge 
interacting with a uniform compensating background 
charge in the periodic cell. Indeed, following the method 
proposed by Ewald^ one replaces the point charge by a 
Gaussian charge distribution, 



1 



8 7 3ti-3/2 



-r 2 /(2 7 ) 2 



(17) 



The second term in Eq. ([9]) then corresponds to the elec- 
trostatic energy of the Gaussian charge distribution in- 
teracting with the background and is evaluated in Fourier 
space. The first term in Eq. ^ compensates for the self- 
interaction of the Gaussian charge. The complete Ewald 
method also gives a third term, which accounts for the 
difference between the point charge and the Gaussian 
charge and which is generally evaluated in real space. 
The absence of this term in Eq. ((9]) accounts for the vari- 
ation of x with 7 in Fig. [TJ However, in the limit 7 — > 0, 
the latter term vanishes and the correction —x/2 pre- 
cisely corresponds to the electrostatic energy of the point 
charge. This connection has been pointed out earlier in 
the case of isolated molecules in large supercells on the 
basis of an intuitive reasoning^ The present derivation 
shows that the same correction term also applies to the 
case of extended systems. 



B. Screened exchange 

We note that the methodology described above does 
not only apply in the presence of a divergence, but might 
also be useful when the interaction potential only shows 
a rapidly varying behavior. Indeed, the direct treatment 
of such a potential in Fourier space is difficult when the 
density of k points cannot easily be increased. For in- 
stance, this applies to screened exchange interactions. 

We here illustrate this point by focusing on the 
screened exchange interaction recently proposed in func- 
tionals developed by Heyd, Scuseria, and Ernzerhof 
(HSE),— but the scheme also applies analogously to other 
forms of screened exchange. In the HSE functional, a 
complementary error function is used to describe the 
short-range exchange interaction: 



V SI = 



erfc(ur) 



(18) 



where u> is a suitable parameter defining the extent of 
the potential. For simplicity, let us again consider the T- 
point approximation. The matrix element of the screened 
exchange operator in a plane-wave basis set is given by 
an expression analogous to Eq. The interaction po- 
tential is given by the Fourier transform of the potential 
defined in Eq. (fl8l). i.e. 



$sr (G) 



1 4tt 



1 



-G 2 /(2w) 2 



(19) 



The G=0 component of this potential <i> sr (G = 0) is not 
divergent, and is equal to ir/(£lui 2 ). However, the lat- 
ter expression cannot be used for any value of u>. In- 
deed, when the screened exchange interaction approaches 
the exact one (viz. in the limit uj — ► 0), <i> sr (G = 0) di- 
verges rather than converging to the correct value given 
in Eq. (fTO]). 

We derive a suitable expression of $ sr (G = 0) by pro- 
ceeding in the same way as in the previous section. The 
correction is again given by Eq. ([7]) , where we now choose 
as auxiliary function: 
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Using the definition given in Eq. @ and taking the limit 
as in Eq. (|10[) . we find the correct expression for the G=0 
component of the screened exchange potential: 



$ sr (G - 0) - X (0) - X 



4w 2 



xH, (21) 



which defines the function x{lo). This leads us to the 
following form for the interaction potential: 

1 4vr 
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(22) 



To illustrate the behavior of this correction as a func- 
tion of u>, we considered a cubic simulation cell with a 
side of 20 bohr. Figure [2a) shows that the function 
x{oj) is nearly indistinguishable from the analytical ex- 
pression ir/(£lui 2 ) when uj > 0.1 bohr -1 . However, for 
lower values of to, the screened potential approaches the 
bare Coulomb interaction and the analytical expression 
gives erroneous results. Instead, the function x(uj) cor- 
rectly reproduces the Coulomb limit for uj — 0. 

For comparison, the parameter uj assumes the value of 
0.106 bohr -1 in the HSE functional^ Hence, for the case 
illustrated in Fig. Ela), the proposed treatment would 
not produce a sizable correction. However, the situation 
changes when smaller simulation cells are used. We show 
in Fig. [2] the dependence of the singularity correction 
on the size of the cubic cell for ui = 0.106 bohr -1 . It 
is seen that the analytical expression deviates from the 
proper correction x f° r cei l sizes smaller than 18 bohr. 
This behavior reflects the fact that at these cell sizes 
the T-point sampling in reciprocal space is too sparse for 
properly treating the spatial varations of the screened 
exchange potential defined by uj = 0.106 bohr -1 . 



III. METHODS 

The semilocal density-functional calculations in this 
work were performed within the generalized gradient ap- 
proximation proposed by Perdew, Burke, and Ernzerhof 
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10 15 20 

L (bohr) 

FIG. 2: (a) The singularity correction for screened exchange 
X [Eq. (|21|l ] as a function of the screening parameter uj for 
a cubic simulation cell of 20 bohr (solid line). At large u 
the singularity correction coincides with the analytical ex- 
pression x(a>) — * 7r/(r2w 2 ) (dashed line), (b) \ vs the side L 
of the cubic simulation cell at fixed to = 0.106 bohr -1 (solid 
line), corresponding to the value set in the Heyd-Scuseria- 
Ernzerhof functional (Ref. [56h . At large L, \ approaches 
ty/(L 3 lo 2 ) (dashed line). The Brillouin zone is sampled at 
the F point. 



(PBE)jl We considered the class of hybrid functionals 
which are obtained by replacing a fraction a of the PBE 
exchange with exact exchange: 57 



£* ybrid = a £™ ct + (l-a)E. 



PBE 



(23) 



In this work, we used the functional defined by a = 0.25, 
which is referred to as PBEOf 5 - 7 " The singularity correc- 
tion then corresponds to a fraction a of the correction 
pertaining to full exact exchange: 
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FIG. 3: (a) Total energies and (b) energy eigenvalue gaps 
of naphthalene calculated using the PBE (left panels) and 
PBEO (right panels) functionals vs supercell size. To allow a 
fair comparison, the same energy scales are used in the left 
and right panels. For PBEO, closed and open symbols indicate 
values obtained with the singularity correction turned on and 
off, respectively. 



carried out with the CPMD code. In this code, we ex- 
tended the available implementation of exact exchange to 
account for the singularity following the scheme outlined 
in Sec. Hi!] 

We also performed electronic structure calculations for 
molecules with the Gaussian03 suite of programs^ We 
used the large cc-pVTZ basis set. This method does not 
give rise to any divergence of the interaction potential 
and allowed us to obtain reliable benchmarks. 



IV. NEUTRAL SYSTEMS 



= a X = 0.25v. 



(24) 



A. Finite systems 



In our electronic structure scheme based on plane- 
waves, only valence electrons are treated explicitly and 
core-valence interactions are described by normconserv- 
ing pseudopotentials i 58 ' 59 Pseudopotentials were gener- 
ated at the PBE level and were also used in the calcu- 
lations based on hybrid functionals. This is expected 
to be a valid approximation for the atoms considered in 
this work, i.e. C, N, O, and Si, in which core-valence 
exchange interactions are weak. We used a kinetic en- 
ergy cutoff of 20 Ry for systems involving only Si atoms, 
but increased the cutoff to 70 Ry when any of the other 
atoms occurred. These cutoffs are sufficiently high to en- 
sure converged total energies and energy eigenvalues. For 
all calculations involving small supercells and dense k- 
point meshes, we used the code PWSCF of the Quantum- 
ESPRESSO package^ in which exact exchange and the 
singularity correction are implemented42, All calculations 
involving large supercells and T-point samplings were 



In this section, we present calculations for total ener- 
gies and energy eigenvalues of small organic molecules. 
For these calculations, we used the electronic structure 
scheme based on plane waves and pseudopotentials with 
large supercells of varying size and T-point sampling. 
We considered the molecules of naphthalene (CioHg) and 
pyridine (C5H5N), which are sufficiently small to allow 
us to investigate the asymptotic behavior for large sim- 
ulation cells. Furthermore, the eigenvalue of the lowest 
unoccupied molecular orbital (LUMO) of these molecules 
is found below the vacuum level, both in PBE and PBEO. 
In Fig. [31 the total energies and the energy eigenvalue 
gaps of naphthalene calculated in PBE and PBEO are 
plotted as a function of the size of the cubic simula- 
tion cell. In hybrid functional calculations, we explic- 
itly compare results obtained with and without the sin- 
gularity correction (3, i.e. obtained by turning on and 
off the G = component of the exact exchange poten- 



TABLE I: Energy eigenvalues of the highest occupied molec- 
ular orbitals (ehomo) and of the lowest unoccupied molecular 
orbitals (elumo), energy eigenvalue band gaps (E g ), and ion- 
ization potentials (IP) of naphthalene and pyridine calculated 
with the PBEO functional through a plane-wave (CPMD, Ref. 
l4ll ) and an all-electron scheme (Gaussian03, Ref. 13(f ). The 
eigenvalues are referred to the vacuum level. 
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tial. One notices that both the total energies and the 
energy eigenvalue gaps calculated in PBEO show a con- 
vergence behavior similar to that achieved in PBE pro- 
vided the singularity is accounted for. Note that for the 
largest considered supercell (side of 60 bohr) the singu- 
larity correction (3 is still quite sizable and amounts to 
about 0.2 eV. A similar behavior is observed for pyri- 
dine (not shown). These results therefore indicate that 
in the case of isolated molecular systems, the singularity 
correction is crucial to achieve well- converged values of 
total energies and energy eigenvalues in hybrid functional 
calculations based on plane- wave basis sets. 

To benchmark our hybrid-functional results, we 
also used an all-electron electronic-structure scheme 
based on localized atomic orbitals and open boundary 
conditions^ Total energies cannot be compared because 
of the different number of electrons that are treated ex- 
plicitly. The energy eigenvalues and energy gaps are com- 
pared with those obtained with the plane-wave scheme in 
Table |TJ The eigenvalues arc referred to the vacuum level 
corresponding to local potential far from the molecule. 
The agreement between the two sets of calculations is 
very good. This agreement further supports the valid- 
ity of the singularity correction derived in Sec. [Til To 
the extent that basis set errors for these molecules are 
small, the differences also provide an estimate of the way 
the different treatments of core electrons affect the en- 
ergy eigenvalues obtained with hybrid functionals. Table 
UJ also contains calculated ionization potentials, of which 
the discussion is deferred to Sec. IV Al 



B. Extended systems 

In this section, we study the effect of the singularity 
correction in hybrid functional calculations on the total 
energies and energy eigenvalues of bulk systems. In par- 
ticular, our purpose is (i) to illustrate the convergence 
of small-cell calculations with fc-point sampling (cf. Refs. 
I461I48T ) and (ii) to study the convergence of T-point cal- 
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FIG. 4: Total energies of (a) Si and (b) Q-quartz Si02 per 
formula unit versus l/jVfc7V at , where is the total number 
of fc-points and N^t the total number of atoms in the super- 
cell. Results obtained in the PBE and the PBEO are reported 
in left and right panels, respectively. For PBEO, closed and 
open symbols indicate values obtained with the singularity 
correction turned on and off, respectively. Arrows show data 
points which were also obtained with F-point sampling. 



culations with supercell size. In the latter case, one issue 
of interest is whether the achieved level of convergence is 
similar to that of corresponding calculations with semilo- 
cal functionals. 

We chose silicon and a-quartz SiC>2 to examine systems 
with different band gaps. To allow a comparison between 
PBE and PBEO calculations of the electronic structure, 
we used the lattice parameters optimized in the PBE also 
in the hybrid functional calculations. 

We first considered the convergence of total energies. 
PBE and PBEO results for silicon and a-quartz are dis- 
played in Fig. |H The data points illustrate how conver- 
gence is achieved with increasing density of fc-point sam- 
pling. For both systems, the inclusion of the singularity 
correction in the hybrid functional calculations leads to 
a faster convergence, very similar to that achieved by the 
semilocal functional. It is apparent that the singularity 
correction is sizable and that its inclusion significantly 
speeds up the convergence. For comparison, we also cal- 
culated total energies employing large supercells and T- 
point sampling. In the respective limits of dense fc-point 
samplings and large supercell size, the two kinds of cal- 
culations give the same converged energies. When the 
energies calculated in the latter scheme are reported in 
Fig. [4] (see arrows) , they are found to correspond to ener- 
gies obtained with primitive cells and dense fc-point sam- 
plings. In particular, for finite supercell size, the degree 
of convergence achieved with hybrid functionals when the 
singularity is treated is comparable to that obtained with 
semilocal functionals. These results highlight the impor- 
tance of including the singularity correction when using 
hybrid functionals with T-point sampling. This is par- 
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FIG. 6: Silicon band gap calculated with the PBE (left panel) 
and PBEO (right panel) functionals for a cubic simulation cell 
containing 8 atoms vs number of fc-points in the [001] di- 
rection, Nkz, with fixed fc-point sampling in the orthogonal 
directions (Nk x =Nk y =2). For PBEO, closed and open sym- 
bols indicate values obtained with the singularity correction 
turned on and off, respectively. 



FIG. 5: Eigenvalues corresponding to the valence (e v ) and 
conduction (e c ) band edges of (a) Si and (b) a-quartz Si02 vs 
l/(NkN at ), where Nk is the total number of fc-points and JV at 
the total number of atoms in the supercell. Results obtained 
in the PBE and the PBEO are reported in left and right panels, 
respectively. For PBEO, closed and open symbols indicate 
values obtained with the singularity correction turned on and 
off, respectively. Arrows show data points which were also 
obtained with T-point sampling. The energies obtained with 
the two functionals are aligned through the local electrostatic 
potential (cf. Ref.[H). 

ticularly important when total energies obtained with 
supercclls of different size are comparatively evaluated, 
such as in the optimization of lattice parameters or in 
constant-pressure molecular dynamics simulations. 

Next, we addressed the convergence of energy band 
edges and band gaps. The band gap is obtained from the 
energy eigenvalues: 

E g = e c - e v (25) 

where e c and e v are the conduction band minimum and 
the valence band maximum, respectively. In Fig. [5J we 
give the band gaps of silicon and SiC>2 as a function of fe- 
point sampling or supercell size, as obtained within both 
the PBE and the PBEO. The results obtained with the 
hybrid functional are given with and without the sin- 
gularity correction. In analogy with the convergence of 
the total energy (Fig. [4|, the convergence of the band 
gap is very slow when omitting this correction. The im- 
provement is more clear for SiC>2, for which the band gap 
converges at a sparser fc-point density. Results obtained 
with large supercells and T-point sampling coincide with 
those obtained with small unit cells and dense fc-point 
meshes. We note that in the case of SiC>2 a converged 
value of the PBEO band gap is already achieved with a 
simulation cell of 72 atoms, well within current compu- 
tationally accessible limits. In the case of Si, the conver- 
gence is slower because of the smaller band gap, but the 
convergence of the hybrid functional result is similar to 



that achieved with the semilocal functional when the sin- 
gularity correction is included. In Fig.[5l the convergence 
of the respective band edges is shown. The singularity 
correction only affects the occupied states by introducing 
a constant downward shift of the energies (cf. Sec. HT|) . 

In many applications concerning surfaces and inter- 
faces, the supercells have an elongated shape to describe 
the transition across the boundary region. For such sys- 
tems, the omission of the singularity correction in a hy- 
brid functional calculation in which a T-point sampling 
is used leads to a peculiar behavior of total energies and 
single particle eigenvalues. To simulate these conditions, 
we considered a cubic 8-atom simulation for bulk silicon 
and increased the fc-point sampling along the [001] direc- 
tion, while keeping the fc-point sampling in the orthogo- 
nal directions constant. This description is equivalent to 
that achieved with an elongated supercell calculation in 
which the Brillouin zone is sampled at the sole T point. 
Figure [5] shows the calculated band gap vs. the number 
of Appoints in the [001] direction. Omitting the singu- 
larity correction leads in this case to a linear increase of 
the band gap. A similar linear increase is also found for 
the total energy (not shown). In neither case, it is there- 
fore possible to obtain a converged value. Including the 
singularity correction reestablishes a converging behavior 
that resembles that found in semilocal density functional 
calculations. 

From the result in Fig. [6l we infer that the singular- 
ity corrections for elongated simulation cells may change 
sign in comparison with cubic unit cells and/or isotropic 
fc-point meshes. To understand this behavior, it is useful 
to recall that, in the case of T-point sampling, the sin- 
gularity correction is proportional to the energy per unit 
cell of a periodically repeated point charge immersed in a 
compensating background. In case of nearly cubic super- 
cells the latter energy is negative, i.e. the attractive inter- 
action of the point charge with the uniform background is 
larger than the repulsive interactions between the point 
charge and its images. When the shape of the cell elon- 
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gates in one or two directions, the repulsive interaction 
with the image charges in the orthogonal directions grows 
because of the reduced screening. For sufficiently elon- 
gated shapes, this repulsive interaction dominates and 
the sign of the point-charge energy switches. 

To summarize the results of this section, we showed 
that the singularity correction is needed to obtain con- 
verged total energies and single-electron eigenvalues 
in hybrid functional calculations with small unit cell 
and dense fc-point samplings, in accord with previous 
studiesj 46 ! 48 Furthermore, we showed that hybrid func- 
tional calculations with large supercells and T-point sam- 
plings also benefit from the singularity correction, yield- 
ing levels of convergence which are similar to those 
achieved with semilocal functionals. In the case of T- 
point samplings, the singularity correction applies to 
both molecular and extended systems in a qualitatively 
similar way. 



V. CHARGED SYSTEMS: LOCALIZED STATES 

In the following, we discuss convergence issues asso- 
ciated to charged systems when using hybrid functional 
schemes based on plane wave basis sets. Charged sys- 
tems occur in several circumstances, such as, for example, 
when studying charged molecules, defects, ions in liquids, 
etc. We here focus on the determination of total energy 
differences between different charge states and the way 
the singularity correction affects the convergence proper- 
ties. We are particularly interested in assessing how the 
convergence properties of hybrid functional calculations 
compare with those of semilocal density functional calcu- 
lations. For simplicity, we consider from now on only sys- 
tems in which the electronic structure is sampled through 
the sole T point. The convergence is therefore studied 
with respect to increasing simulation cell, or equivalently 
with respect to decreasing singularity correction. Gener- 
alization to convergence with /c-point samplings is trivial. 
In this section, we deal with atomically localized states, 
either in finite systems or as defect states in solids. In- 
finitely delocalized states are discussed in Sec. I VII In 
this work, we do not consider states showing intermedi- 
ate degrees of localization. For a discussion of the latter, 
we defer the reader to Refs. l6ll 



A. Finite systems 

To illustrate the convergence of total energy differences 
between different charge states, we considered the cal- 
culation of the ionization potential (IP) of an isolated 
molecule: 

IP = E N -i — Eff, (26) 

where Em is the total energy of the neutral molecule and 
En—i the total energy of the same molecule in which 
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FIG. 7: Ionization potential (IP) of naphthalene calculated 
with the PBE (left panel) and PBEO (right panel) functionals 
for cubic simulation cells vs singularity correction f3, which 
scales like the inverse of the simulation cell size. For PBEO, 
closed and open symbols indicate values obtained with the 
singularity correction turned on and off, respectively. 

one electron has been removed. The total energy cal- 
culations correspond to an isolated molecule placed in 
a large supercell and subject to periodic boundary con- 
ditions. This technical constraint introduces a spurious 
interaction between the localized charge and the neutral- 
izing background charge, which needs to be considered in 
order to speed up the convergence^ This correction de- 
pends on the size of the simulation cell, but applies indif- 
ferently to both hybrid functionals and semilocal density 
functionals. In our calculations, the dominant correction 
corresponding to the charge monopole has systematically 
been included.— 

In Fig. [71 we give the ionization potentials of naph- 
thalene calculated within both the PBE and the PBEO 
for different supercell sizes. The PBEO results are ob- 
tained by both including and dismissing the singularity 
correction. The results are plotted as a function of the 
singularity correction, which scales like the inverse of the 
simulation cell size. For the simulation cells considered 
(cubic cells with sides ranging between 20 and 60 bohr), 
the PBE results are already close to the converged values 
obtained by linear extrapolation. The same consideration 
applies for the PBEO results which include the singular- 
ity correction. However, in absence of singularity cor- 
rection, the error with respect to the converged values is 
significantly larger. The converged PBEO values for the 
ionization potentials of naphthalene and pyridine are re- 
ported in Table [H where they are compared with results 
obtained with an all-electron scheme based on localized 
orbitals. The values calculated within the two schemes 
differ by less than 0.05 eV. 

It is of interest to extend our comparative study be- 
tween PBE and PBEO calculations to the approximate 
scheme based on Slater's transition stated According to 
the integral form of Janak's theorem,— the total-energy 
difference of Eq. (|2"6")) is given by 

E N _ 1 -E N = - [ e n (f)df, (27) 
Jo 

where e n (f) describes the eigenvalue of the highest oc- 
cupied eigenstate n as it varies with its fractional occu- 
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pation /. Using the trapezoidal rule for the integral, we 
obtain Slater's approximation: 

E N -! -E N & -e„(l/2) « - [e„(0) + e„(1)] /2, (28) 

where we further approximated the eigenvalue at half- 
filling by an average at integer occupations. For the 
molecules investigated here, this approximation is found 
to give accurate results for both the semilocal and the 
hybrid functional (not shown). Note that, in the latter 
case, the singularity corrections of energies and eigenval- 
ues are compatible with the relation in Eq. ([28]) . Indeed, 
this approximation equally holds for corrected energies 
and eigenvalues as for uncorrected ones. Similar consid- 
erations apply to electron removal energies. 



B. Defects in solids 

The discussion pertaining to total energy differences 
of finite systems (Sec. IV A[) applies with minor modifi- 
cations to the study of localized defect states in solids. 
In this case, the relevant physical quantities, the charge 
transition levels, are also expressed as total energy dif- 
ferences. Defect formation energies are first determined 
for varying electron chemical potential //jiS 
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(29) 

where E q ot is the total energy of the defect system car- 
rying a charge q, S^ot the total energy of the unper- 
turbed system, n a the number of extra atoms of species 
a needed to create the defect, and r\ a the correspond- 
ing atomic chemical potential. The chemical potential 
/i is referred to the valence band maximum e v . AV is 
a correction which is applied in order to align the local 
potential far away from the neutral defect to that of the 
unperturbed bulkji& The correction E q OTT describes the 
spurious interaction of the added charge with the com- 
pensating background charge. 62 As commented in Sec. 
Ill A[ the leading term of this correction pertaining to the 
monopole can be expressed as 



E q 



q 2 x 

2e ' 



(30) 



where \ is defined in Eq. ([TO")) and e is the dielectric con- 
stant of the unperturbed bulk system. In our calcula- 
tions, this correction is applied systematically in both 
PBE and PBEO calculations. 

Charge transition levels correspond to specific values 
of the electron chemical potential for which two charge 
states have equal formation energies. For example, the 
charge transition level n q / q ' between two charge states q 

and q' is defined by the condition E q — E q and is thus 
given by the following expression: 
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FIG. 8: Vertical charge transition levels /i + /o associated to 
the hydrogen bridge defect in /3-crystobalite as obtained with 
the PBE (left panel) and the PBEO (right panel) functionals 
vs number of atoms JV at in the simulation cell. For PBEO, 
closed and open symbols indicate values obtained with the 
singularity correction turned on and off, respectively. The 
indicated valence and conduction band edges correspond to 
the converged energy eigenvalues. The energies obtained with 
the two functionals are aligned through the local electrostatic 
potential (cf. Ref.|H). 



In this expression the dependence on the atomic chemi- 
cal potentials r] drops out and the defect charge transi- 
tion level is basically determined by a total energy differ- 
ence between different charge states of the defect. In this 
sense, these quantities are counterparts of the ionization 
potentials and electron affinities of isolated atomic and 
molecular systems. The charge transition levels of local- 
ized defect states can also be obtained in a very accurate 
way through the energy eigenvalues by the application of 
Janak's theorem both in PBE and in PBEO [Eq. (|2^)) ]. 65 

We illustrate the convergence behavior of charge tran- 
sition levels for the hydrogen bridge defect (hydrogen 
substitutional to oxygen) in /3-crystobalite SiC-2. In par- 
ticular, we focused on the vertical charge transition level 
/i + /o located at approximately mid gapi 66 ' 67 We consid- 
ered three different supercells in which one lateral size 
was varied while the other dimensions were kept fixed, 
i.e. we used 2x2x2, 2x2x3, and 2x2x4 unit cells. Such 
elongated cells might for instance occur when slab models 
are adopted. Charge transition levels calculated in the 
PBE and in the PBEO are shown in Fig. [5] for increasing 
number of atoms in the simulation cell. The results ob- 
tained with the two functionals are aligned through the 
local electrostatic potential, as suggested by the study in 
Ref.[29|. It clearly appears that the PBE and PBEO defect 
levels show a similar convergence behavior, provided the 
singularity correction is included in the PBEO calcula- 
tion. The singularity correction is crucial to achieve this 
level of convergence. Without the singularity correction, 
the charge transition levels are clearly not converged for 
the range of simulation cells considered. It is clearly seen 
that the convergence behavior of charge transition levels 
of defects is analogous to that of energy transitions in 
finite systems fSec. IV A"|) . 
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VI. CHARGED SYSTEMS: 
STATES 



DELOCALIZED 



PBE 
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The case in which the extra charge is carried by ex- 
tended delocalized states requires special attention. Let 
us assume an infinite solid with a energy eigenvalue spec- 
trum £fc. The energy cost of adding one electron to the 
previously unoccupied state n is simply given by 



En+i — En = e 7 



(32) 



where we used Janak's theorem^ as formulated in Eq. 
(|27| and the fact that energy eigenvalues in infinite solids 
do not depend on the occupation of the state. Identical 
considerations also apply for electron removal energies. 
On this basis, it is possible to express the valence and 
conduction band edges in terms of total energy differ- 
ences between systems of different charge: 



e c — En+i — En, 
e v = E N — En-i- 



(33) 
(34) 



Consequently, the energy band gap can also be obtained 
as 
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(35) 



In our notation, the tilde sign signifies that the concerned 
quantity is obtained from a total energy difference, as op- 
posed to a direct derivation from the spectrum of energy 
eigenvalues. 

In practical calculations, the supercells always have 
finite size, by which the band edges determined by 
total-energy difference differ from the actual energy 
eigenvalues^ 
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FIG. 9: Band edges of a-quartz calculated with the PBE 
(left) and PBEO (right) functionals as total energy differences 
vs singularity correction /3 that characterizes the simulation 
cell, the limit of infinite simulation cell being achieved for 
vanishing /3. For PBEO, closed and open symbols indicate 
values obtained with the singularity correction turned on and 
off, respectively. The energies obtained with the functionals 
are aligned as in Fig. [5] 
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FIG. 10: Difference between band edges of a-quartz evaluated 
through total energy differences (e c , e v ) and as energy eigen- 
values (e c , £ v ) in the PBE (left panel), (b) in the PBEO with 
singularity turned off (middle panel), and (c) in the PBEO 
with the singularity correction turned on (right panel). The 
results are plotted against the singularity correction /3. In the 
right panel, the dashed lines correspond to Ae c = — /3/2 and 
Ae v = +,3/2. 



where f2 is the volume of the simulation cell. Similarly, 
we define 



E g = E s + AE g , 



(40) 



where AE g — Ae c — Ae v . 

We first focus on results obtained with the semilocal 
density functional. For illustration, we considered a- 
quartz SiC>2 and used simulation cells of varying size. In 
Fig. [51 we report the band edges calculated using the total 
energy differences given in Eqs. (|3"3")) and (f3"4"|) . and com- 
pare them with the converged energy eigenvalues. For 
the simulation cells considered here, the band edges ob- 
tained in the two different ways are essentially identical, 
i.e. Ae c w and Ae v w 0. 



Figure [H also shows corresponding results obtained 
with hybrid functionals. In contrast to the behavior 
found for semilocal functionals, we now encounter a much 
slower convergence. In other words, in hybrid functional 
calculations, Ae c and Ae v are significantly larger than in 
semilocal density functional calculations performed with 
the same simulation cells. Therefore, we infer that this 
behavior should be ascribed to the singular behavior of 
the exchange interaction. 

In order to understand this behavior, we reason as fol- 
lows. The relations given by Eqs. (|36)) and (j37|) follow 
from the Janak theorem and apply to any analytical func- 
tional. Therefore, they also apply to a case functional 
which is specific to a given simulation cell and a given 
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basis set. We further define the case functional by set- 
ting the G = component of the exchange potential to 
zero, which corresponds to the label "uncorr" introduced 
earlier. For this case functional, Ae c and Ae v are ex- 
pected to behave in a similar way as for a semilocal func- 
tional. Hence, for simulation cells large enough to yield 
vanishing Ae c and Ae v with the semilocal functional, we 
similarly expect vanishing Ae" ncorr and Ae" ncorr . This 
implies 



puncorr puncorr 

jriuncorr puncorr 
~ &N-1 ■ 



(41) 
(42) 



As shown in the first two panels in Fig. [TQ1 this behavior 
is numerically confirmed for our case study of a-quartz. 

Using the relationship between corrected and uncor- 
rected quantities determined in Eqs. (Tl3]) and (|16p . we 
then derive: 
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In other words, hybrid functional calculations in finite 
simulation cells give 



Ae c 
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(45) 
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This result is graphically illustrated for a-quartz in the 
third-column panels of Fig. [TUJ For the band gap, this 
implies: 



(47) 



We note that the remaining dependence of both e c and 
e v on (3 does not result from the integration of the sin- 
gularity, since the results in the third panel in Fig. [TU] 
already refer to "corrected" results. The remaining dif- 
ferences rather originate from the exchange selfinterac- 
tion associated to the extra charge, which vanishes slowly 
with increasing simulation cell. Hence, this behavior is 
not specific to the use of plane- wave basis sets and should 
also manifest in implementations based on other basis 
functions. 

Our findings concerning the energy eigenvalues are 
schematically illustrated in Fig. QT] The scheme refers 
to a case in which all quantities would already be con- 
verged with simulation cell size, if it were not for the 
occurrence of exact exchange. Figure [TlT a) refers to the 
corrected energy eigenvalues and corresponds to the con- 
verged result. Figure [TlTb') shows the eigenvalues in the 
"uncorrected" case. The singularity correction shifts the 
occupied states downwards by (3, leaving the unoccupied 
ones unaffected. Figure [TlTc) refers to the determination 
of band edges through the use of total-energy differences. 
It is seen that the band edges obtained in this way still 
differ from the converged levels, despite the use of "cor- 
rected" quantities. The valence band edge is overesti- 
mated by (3/2, whereas the conduction band edge is un- 
derestimated by (3/2. Consequently, the band gap is un- 
derestimated by (3 through this approach. Figure fTUd) 
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FIG. 11: Schematic representation of the band edges deter- 
mined with hybrid functionals in the presence of a limited k 
point sampling: as energy eigenvalues (a) with and (b) with- 
out accounting for the singularity correction and through to- 
tal energy differences (c) with and (d) without singularity 
correction. The result in (a) shows the fastest convergence 
with simulation cell. 



also refers to band edges obtained through total energy 
differences but without including the singularity correc- 
tion. One obtains the same result as in Fig. [TlT b). illus- 
trating thereby that the integration of the singularity in 
the exchange term is responsible for the slow convergence 
of the band edges calculated by total-energy differences. 



VII. CONCLUSIONS 

In this work, we investigated the use of hybrid func- 
tionals in plane- wave implementations, in comparison 
with semilocal density functionals. The main objective 
consisted in determining whether a hybrid functional cal- 
culation requires a fc-point sampling of increased density 
in order to properly integrate the singularity appearing 
in the exact nonlocal exchange energy. This issue is of 
particular importance when dealing with large simulation 
cells, where an excessive increase of the fc-point sampling 
would make the calculation computationally prohibitive. 
Typical applications include surface, interface, and defect 
calculations, but also molecular dynamics simulations. 

To treat the divergence, we adopted a formulation 
which consists in transforming the integrand into a reg- 
ular function through the use of an auxiliary function 
that can be integrated analytically^ Through the use 
of an appropriate auxiliary function— this formulation 
can trivially be extended to calculations with large sim- 
ulation cells and low-density fc-point samplings. In the 
case of T-point sampling, the sampling of reciprocal space 
is achieved through the reciprocal lattice vectors, which 
densify as the simulation cell grows. We further used a 
formulation which recasts the treatment of the divergence 
in the form of singularity correction terms.— These terms 
intervene in the total energy and in the energy eigenval- 
ues of the occupied states. 

In the present investigation, we found it convenient 
to distinguish finite and extended systems, localized and 
delocalized states, and neutral and charged calculations. 
The general conclusion is that the same /c-point sam- 
plings used in semilocal density functional calculations 
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yield a comparable level of convergence in hybrid func- 
tional calculations, provided the singularity corrections 
are accounted for. However, our study highlights a few 
points that deserve special attention. 

The first point concerns the treatment of screened ex- 
change. While screened exchange does not show any 
singularity, it is nevertheless recommended to adopt the 
proposed scheme also in this case in order to achieve the 
same convergence properties as achieved with semilocal 
density functionals. 

The second noteworthy point concerns applications in 
which the sampling in reciprocal space around the di- 
vergence is anisotropic. This is for instance the case for 
calculations with elongated simulation cells and T-point 
samplings, as often occurs in the study of surfaces and 
interfaces. In such cases, the singularity correction terms 
are critical for achieving not only well converged proper- 
ties but also a qualitatively correct behavior. 

The final point that deserves attention and which is 
unusual with respect to ordinary semilocal density func- 
tional calculations concerns the determination of band 
edges and band gaps of extended systems through total 
energy differences. Our study shows that band edges de- 
termined in this way converge slower in hybrid functional 



calculations because of the exchange selfinteraction asso- 
ciated to the extra charge. This convergence problem 
arises for delocalized states but does not occur for local- 
ized states of point defects or of finite molecular systems. 

In conclusion, the correct treatment of the singular- 
ity can be achieved without requiring any significant 
computational overhead. This opens the way for us- 
ing hybrid functionals in very much the same way as 
ordinary semilocal functionals. To date, the scheme de- 
scribed in this work has already led to several successful 
applications including studies of amorphous systems,- 8 
defect8^S^2mS2 and interfaces.2S^2JiZLE2^ 
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